function varR2 = model_varR2(w,tspan,dt)
% varR2 = model_varR2(w,tspan)
% R^2 = \sum_{i=1}^{MODE}\alpha_i^2+\beta_i^2
% var(R^2) = sum_{i=1}^{MODE}var(\alpha_i^2)
%            + sum_{i=1}^{MODE}var(\beta_i^2)
% 
varR2 = zeros(length(tspan),1);
varAlpha2 = zeros(MODE,1);
varBeta2  = zeros(MODE,1);

% Calculate variance, which is the diagonal elements of the covariance matrix
for i =1:MODE
    temp1 = exp(-2*nu*i^2*dt);
    temp2 = exp(-4*nu*i^2*dt);
    temp3 = sigma^2/(nu*i^2);
    temp4 = temp3*(2*temp4-2*temp3);
    temp5 = 0.5*temp3^2*(temp2-2*temp1+1);
    varAlpha2(i) = -(alpha0(i)^2)*temp4+temp5;  
    varBeta2(i) = -(beta0(i)^2)*temp4+temp5;
end

varR2 = sum(varAlpha2+varBeta2);